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I.  Introduction  and  Existing  State-of-the-Art 

Vehicles  traveling  through  an  atmosphere  at  hypersonic  speeds  generate  strong  shock  waves.  These  shock 
waves  generate  extremely  high  gas  temperatures  within  the  shock  layer  between  the  shock  itself  and  the 
surface  of  the  vehicle.  As  the  gas  reaches  such  high  temperatures,  its  vibrational  degrees  of  freedom 
become  excited  and  diatomic  and  polyatomic  molecules  dissociate  into  reactive  atomic  species.  At  high 
altitudes  characteristic  of  hypersonic  flight  where  the  free-stream  density  is  low,  these  reactive  species 
diffuse  through  the  boundary  layer  and  may  chemically  react  with  the  vehicle’s  thermal  protection  system 
(TPS).  Many  TPS  materials  act  as  a  catalyst  for  the  heterogeneous  recombination  of  dissociated  species 
back  into  molecules.  Such  exothermic  surface  reactions  transfer  additional  energy  to  the  vehicle  surface 
and  thus  contribute  to  the  heat  flux  and  overall  heat  load  that  the  TPS  must  withstand.  For  example, 
studies  have  shown  that  surface  catalysis  can  contribute  up  to  30%  of  the  total  heat  load  for  Earth  reentry 
[1]. 

Typically  at  the  beginning  of  the  design  phase  of  a  TPS,  100%  catalytic  efficiency  is  assumed  at  the 
vehicle  surface.  This  provides  an  upper-limit,  conservative  estimation  for  the  heating  rates.  For  100% 
catalytic  efficiency,  all  dissociated  atoms  impacting  the  surface  are  assumed  to  leave  the  surface  as 
recombined  molecules.  Later  in  the  design  cycle  the  boundary  conditions  may  be  relaxed  to  a 
parameterized  catalytic  efficiency  (using  a  recombination  rate  y,  that  is  a  function  of  temperature  and 
pressure)  if  such  data  exists  for  the  specific  gas  and  solid  conditions  involved.  As  outlined  in  a  recent 
review  article  on  the  status  of  aero  thermal  modeling  [2],  such  data  is  very  limited  and  contains  large 
uncertainties.  While  lower  and  upper  bounds  on  heating  rates  are  certainly  useful,  depending  on  the  flight 
conditions  and  atmosphere,  these  predictions  can  differ  significantly.  For  example,  the  stagnation  point 
heat  flux  for  Mars  entry  varies  by  a  factor  of  three  between  the  assumptions  of  a  highly-  and  weakly- 
catalytic  wall  [3].  An  example  of  experimental  data  obtained  on  RCG  coated  Space  Shuttle  tiles  [4]  is 
shown  below  in  Fig.  1.  Here  the  atom  recombination  coefficient  y  (for  both  oxygen  and  nitrogen  atoms)  is 
curve-fit  to  various  functions  of  temperature. 


Nitrogen  Oxygen 


l/T  [K1]  1/T  [K1] 


Figure  1  -  Experimentally  determined  recombination  rate  (y)  as  a  function  of  temperature  for  RCG  coated  tiles 
(taken  from  [4]  “Surface  Catalysis  and  Characterization  of  Proposed  Candidate  TPS  for  Access-to-Space  Vehicles”, 
Stewart  1997). 

Currently  there  is  no  fundamental  explanation  for  these  complex  trends  with  temperature.  It  is  possible 
the  large  changes  in  slope  are  due  to  different  chemical  mechanisms  being  “activated”  at  certain 
temperature  thresholds,  but  these  mechanisms  are  unknown.  Inherently,  the  parameter  y  represents  the  net 
effect  of  many  elementary  chemical  mechanisms  that  collectively  lead  to  recombination.  As  a  result,  it  is 
difficult  to  relate  variations  in  y  directly  to  fundamental  surface  chemistry.  This  inhibits  a  true 
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understanding  of  the  phenomenon  and  limits  the  transferability  of  the  model.  Furthermore,  the  high 
temperature  data,  obtained  in  arc -jet  tests,  are  inferred  from  heat  flux  measurements.  This  requires  an 
accurate  model  of  the  arc -jet  flow  field  and  predicted  convective  heat  flux.  When  compared  with  the 
measured  heat  flux,  the  discrepancy  is  attributed  to  recombination,  where  full  energy  accommodation  is 
also  assumed.  In  fact,  the  arc -jet  flow  is  poorly  characterized  and  the  assumption  of  full  energy 
accommodation  is  not  necessarily  valid.  Thus,  extrapolation  of  these  inferred  y  values  to  flight  conditions 
is  not  rigorous. 


II.  Prior  Research  Efforts 

II.  1  Finite-Rate  Gas-Surface  Models 

A  more  fundamental  model  for  gas-surface 
chemistry  is  analogous  to  the  gas-phase  finite- 
rate  chemistry  models  that  have  been 
successfully  used  for  many  years.  However, 
in  practice,  the  parameterization  of  such 
finite-rate  gas-surface  interaction  models  is 
extremely  challenging  [1,5-14].  For  gas- 
surface  systems,  the  chemical  reactions  now 
involve  gas-phase  species  as  well  as  adsorbed 
species,  and  reactions  should  be  made 
specific  to  various  “chemical  sites”  on  the 
surface.  As  an  example,  a  test  model  for  air  is 
listed  in  Table  1.  The  model  contains  air 
species  N2,  02,  N,  O,  NO,  as  well  as  adsorbed 
atomic  species  denoted  as  A(s),  metastable 
(diffusing)  adsorbed  species  A(s)m,  and 
finally  open  surface  sites  [s].  The  reaction 
mechanisms  include  adsorption/desorption 
reactions  as  well  as  Eley-Rideal  (E-R)  and 
Langmuir-Hinshelwood  (L-H)  recombination 
reactions.  E-R  reactions  occur  when  a  gas- 
phase  atom  strikes  an  adsorbed  atom  and 
recombines  into  a  molecule  that  leaves  the 
surface.  L-H  reactions  occur  when  two 
adsorbed  atoms  diffuse  into  close  proximity 
and  recombine  to  form  a  molecule  that  leaves 
the  surface.  Reverse  reactions  of  these  processes  are  also  included. 
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Table  1  -  Gas-surface  reactions  for  a  generic  air  surface 
catalysis  model. 


The  net  flux  mediated  by  the  surface  depends 
on  the  dynamics  of  each  species  on  the  active 
surface  sites  according  to  the  mechanisms 
listed  and  on  the  rate  constants  associated 
with  each  mechanism.  For  example,  the 
specific  form  of  the  E-R  rate  constant  is 
shown  in  Table  2  with  an  initial  qualitative 
parameterization  [5].  Currently  there  is  no 
way  of  experimentally  determining  all  the  required  parameters  and  rate  coefficients.  However,  these 
parameters  now  have  a  direct  physical  link  to  fundamental  chemistry.  Specifically,  the  coefficients 
correspond  to  activation  energies  (£’),  steric  factors  (y),  surface  site  concentrations  (Ototai),  and  kinetic  flux 
quantities  (v/4).  In  principle,  the  reaction  mechanisms  and  parameters  required  for  their  rates  could  be 


Eley  -  Rideal 
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Table  2  -  Example  of  a  test 
parameterization  for  the 
forward  reaction  rates. 


3 


determined  by  computational  chemistry  studies.  Indeed,  this  approach  has  been  successful  for  many  gas- 
phase  finite-rate  chemistry  models  currently  in  use  and  our  research  attempts  to  accomplish  this  for  gas- 
surface  finite-rate  chemistry. 

II.2  Computational  Chemistry  Studies  of  Unphysical  Crystalline  Surfaces 

A  substantial  amount  of  previous  research  had  been  conducted  in  this  area  by  two  European  research 
groups;  one  located  in  Bari,  Italy  [15-16]  and  the  other  in  Barcelona,  Spain  [17-21].  In  fact,  the  group 
from  Spain  recently  published  an  article  (in  2011)  where  a  finite-rate  catalytic  model  is  developed  for  air- 
silica  interactions  using  computational  chemistry  simulations  [21].  We  have  performed  a  thorough  review 
of  this  prior  work  since  it  is  so  closely  related  to  the  main  objectives  of  this  AFOSR  grant.  Unfortunately, 
during  the  early  stages  of  our  research  for  this  grant,  we  quickly  realized  that  these  prior 
methodologies  were  completely  inadequate  for  real  materials  used  in  hypersonic  applications,  and 
we  found  the  conclusions  of  many  of  these  studies  to  be  invalid.  Since  this  AFOSR  research  report  will 
present  a  new  oxygen-silica  finite  rate  catalytic  model  that  is  substantially  different  from  this  previously 
published  work,  a  critical  review  of  this  prior  work  is  warranted.  A  brief  review  follows  that  clearly 
shows  a  large  disconnect  between  the  computational  chemistry  community  and  the 
aerothermodynamics  community  for  these  problems  and  highlights  some  of  the  critical  issues 
addressed  by  our  AFOSR  grant  research. 

Unphysical  Bulk  Material  (Prior  Studies): 

Both  research  groups  used  computational  chemistry  techniques  to  investigate  oxygen  interactions  with  a 
specific  crystalline  polymorph  of  Si02  (called  (3-cristobalite).  Computer  images  of  this  crystal  lattice  are 
shown  in  Fig.  2.  The  choice  of  (3-cristobalite  is  motivated  by  experimental  studies  from  Balat-Pichelin  et 
al.  [15,22,23]  where  silicon-carbide  (SiC)  surfaces  were  exposed  to  high  temperature  air-plasmas.  During 
such  exposure,  an  oxide  layer  was  quickly  formed  and  atomic  oxygen  loss-rates  (recombination 
coefficients)  were  then  determined  next  to  the  oxide  layer  surface.  While  the  experimental  measurements 
are  of  high  quality  and  repeatability,  the  conclusion  that  the  measured  loss  rates  correspond  to  |3- 
cristobalite  was  based  on  the  fact  that  for  polymorph  diagrams  of  Si02,  (3-cristobalite  is  most  stable 
within  the  experimental  surface  temperature  range.  However,  it  is  almost  impossible  to  imagine  that  a 
crystalline  oxide  layer  was  formed  by  exposing  SiC  to  an  air  plasma.  Indeed,  production  of  crystalline 
materials  (such  as  quartz,  another  stable  Si02  polymorph)  requires  careful  and  precise  manufacturing 
techniques  to  slowly  grow  the  crystal  without  introducing  defects.  Thus,  these  prior  studies  are 
computationally  investigating  a  different  material  than  used  in  the  experiments  they  compare  with. 
Rather,  the  oxide  layer  produced  in  the  experiments  is  most  likely  amorphous  Si02. 

Unphysical  Surface  Structure  (Prior  Studies): 

Furthermore,  in  order  to  simulate  interactions  with  a  surface ,  these  prior  studies  simply  “cleaved”  the 
bulk  crystal  along  the  (001)  plane  (as  seen  in  Fig.  2a)  and  then  proceeded  to  study  gas-phase  collisions 
with  this  surface.  However,  such  cleaved  surfaces  have  highly  unphysical  dangling  bonds  (see  the  top 
layer  of  Si  atoms  in  Fig.  2a).  For  example,  these  studies  conclude  (and  we  have  confirmed)  that  any  gas- 
phase  oxygen  that  hits  these  dangling  bonds  will  immediately  adsorb  (as  seen  in  Fig  2b)  with  such  a 
strong  bonding  energy  that  no  future  catalytic  reaction  could  possibly  remove  this  atom  from  the  surface 
to  form  a  molecule. 
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(a)  Silicon  Terminated  (001)  0-  (b)  Oxygen  Terminated  (001)  0- 

Cristobalite  Surface  (side)  Cristobalite  Surface  (side).  Topmost 

layer  of  oxygen  highlighted  in  blue. 

Figure  2  -  Cleaved  (001)  (3-cristobalite  surfaces  used  in  prior  computational  chemistry  studies  of  surface  catalysis 
under  conditions  relevant  to  hypersonic  flight  (Si  -  yellow,  O  -  red,  O  originating  from  the  gas  are  blue). 

In  order  to  study  recombination,  these  previous  studies  instead  placed  adsorbed  oxygen  atoms  at  different 
(and  less  energetically  bound)  locations  and  then  studied  collisions  with  incoming  gas-phase  oxygen 
atoms.  This  now  led  to  a  certain  rate  of  recombination  reactions,  which  of  course  depended  upon  where 
the  adsorbed  oxygen  was  placed.  By  following  this  invalid  modeling  approach,  remarkable  agreement 
with  measurements  from  air  plasma  experiments  was  reported  [15,16].  In  contrast,  our  computational 
chemistry  simulations  predict  that  such  cleaved  surfaces  reconstruct  into  the  lower-energy,  stable  surface 
configurations  seen  in  Fig.  3.  Such  surface  reconstructions  are  supported  by  a  number  of  experimental 
studies  that  will  be  detailed  later  in  this  report.  Thus,  not  only  is  the  bulk  material  ((3-cristobalite) 
unphysical  but,  perhaps  more  importantly,  the  surfaces  that  gas-phase  atoms  interact  with  in  these 
previous  studies  are  not  only  physically  unrealistic,  but  they  are  computationally  unstable. 


(a)  Silicon  Terminated  (001)  0-  (b)  Oxygen  Terminated  (001)  0- 

Cristobalite  Surface  (side)  Cristobalite  Surface  (side).  Topmost 

layer  of  oxygen  are  highlighted  in  blue. 


Figure  3  -  Surface  reconstmctions  of  (001)  (3-cristobalite  (Si  -  yellow,  O  -  red,  O  originating  from  the  gas  are  blue). 
Arbitrary  Surface  Coverage  (Prior  Studies): 

Once  probabilities  of  adsorption,  desorption,  and  recombination  are  determined  in  function  of  surface 
temperature  and  gas-phase  impact  energy,  one  can  develop  a  reaction  rate  model  parameterized  in 
function  of  the  gas-surface  temperature  (T).  If  these  rate  equations  are  solved  for  a  given  partially- 
dissociated  gas  state  (species  concentrations,  temperature,  and  pressure),  the  resulting  concentrations  of 
adsorbed  species  [A(s)]  can  be  calculated  from  the  reactions  in  Table  1.  These  adsorbed  concentrations 
then  correspond  to  a  specific  surface  coverage  corresponding  to  these  specific  gas  conditions. 
Determination  of  the  surface  coverage  is  extremely  important,  as  it  dictates  the  actual  surface  that 
gas-phase  molecules  interact  with.  In  the  previously  published  finite  rate  model  [21],  the  model 
predicted  a  steady-state  coverage  saturated  with  02  molecules  for  which  no  further  recombination  is 
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possible.  However,  a  surface  coverage  part-way  between  zero  coverage  and  steady  state  was  selected  and 
recombination  efficiencies  for  this  arbitrary  coverage  were  reported,  which  also  happened  to  match 
experimentally  determined  values  very  well.  There  is  no  basis  for  assuming  such  an  arbitrary  surface 
coverage.  In  fact,  a  steady  state  surface  coverage  is  established  over  very  small  time  scales  while  catalytic 
reactions  continue  on  the  surface  for  macroscopic  time  scales.  Thus,  the  surface  coverage  assumed  in 
prior  studies  was  not  representative  of  the  experiments  they  compared  with,  and  was  even  in 
contradiction  with  the  computational  chemistry  simulations  of  the  same  studies. 

Improper  Application  of  Detailed  Balance  and  Unphysical  Backwards  Rates  (Prior  Studies): 

Finally,  as  will  be  shown  from  our  research,  when  forming  a  finite-rate  chemistry  model,  detailed 
balance  must  properly  be  accounted  for  such  that  reverse  reaction  rates  and  equilibrium  states 
remain  physically  realistic.  Specifically,  in  the  previously  published  rate  model  [21],  the  reverse  rate  of 
E-R  recombination  (dissociative  adsorption  of  a  gas-phase  molecule)  is  very  high  and  roughly  equal  to 
that  of  atomic  oxygen  adsorption.  However,  the  02  surface  impact  simulations  presented  in  the  same 
publication  do  not  find  any  dissociative  adsorption  reactions  despite  finding  an  atomic  adsorption 
probability  close  to  one.  Furthermore,  the  reported  desorption  rate  coefficient  is  higher  than  the 
theoretical  limit  (set  by  the  vibrational  frequency  of  02  molecules)  by  two  orders  of  magnitude.  Thus,  this 
previously  published  finite-rate  catalytic  model  for  oxygen-silica  interactions  contains  backwards 
rates  that  are  both  unphysical  and  in  contradiction  with  the  computational  chemistry  simulations 
themselves. 

These  prior  studies  [15-21]  have  been  published  in  high  quality  computational  chemistry  and  surface 
science  journals.  The  above  critical  review  clearly  demonstrates  the  lack  of  rigorous  collaboration  and 
understanding  between  computational  chemistry  researchers  and  aerothermodynamics  researchers. 
Specifically,  the  conclusions  of  our  research  under  this  AFOSR  grant  are  that  each  of  the  above  aspects, 
bulk  material,  surface  structure,  surface  coverage,  and  detailed  balance  are  crucial  to  developing  a 
physically  realistic  and  consistent  finite  rate  model  for  oxygen-silica  catalytic  reactions. 

II. 3  Computational  Chemistry  Studies  of  Amorphous  SiO?  Surfaces 

Co-PI,  Dr.  Cozmuta,  published  an  AIAA  conference  paper  in 
2007  entitled,  “Molecular  mechanisms  of  gas  surface 
interactions  in  hypersonic  flow”  [24].  In  this  paper,  Cozmuta 
presented  preliminary  results  for  computational  chemistry 
simulations  of  dissociated  air  interacting  with  an  amorphous 
Si02  surface.  The  simulations  began  with  a  vacant  amorphous 
Si02  surface  exposed  to  partially-dissociated  high  temperature 
air.  Recombination  reactions  were  monitored  only  after  the 
surface  coverage  reached  a  steady  state.  An  image  of  such  a 
simulation  is  shown  in  Fig.  4.  This  type  of  simulation  is  much 
more  physically  realistic  than  those  performed  in  the  studies 
reviewed  above,  and  co-PI  Cozmuta’s  work  provided  the 
original  framework  for  our  AFOSR  grant  research. 

Specifically,  Cozmuta  chose  to  use  a  reactive  force-field  to 

model  inter-atomic  forces,  called  ReaxFF.  The  unique  ReaxFF 
formulation  enables  large  scale  Molecular  Dynamics  (MD) 
simulation  while  accurately  modeling  chemical  reactions. 

ReaxFF  will  be  described  in  detail  in  an  upcoming  section  of 
this  report. 


Figure  4  -  ReaxFF  MD  simulations  of 
dissociated  air  interacting  with  amorphous 
Si02. 
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III.  AFOSR-Grant  Research  Accomplishments 

Overview  of  main  accomplishments 

The  accomplishments  are  subdivided  into  the  following  four  main  subsections: 


III.  1  Implementation  of  a  Finite-Rate-Catalytic  Boundary  Condition  in  the  US3D  CFD  Code 

a  new  finite-rate  catalytic  (FRC)  wall  boundary  condition  has  been  implemented  in  the  US3D 
CFD  developed  and  maintained  at  the  University  of  Minnesota 

this  FRC  boundary  condition  solves  a  set  of  user-specified  chemical  reactions  at  all  wall 
elements  in  a  general  3D  unstructured  CFD  simulation 

the  implementation  has  been  extensively  validated  with  an  independent  implementation  of  the 
same  technique  in  the  NASA  DPLR  code  through  code-to-code  comparisons 
the  behavior  of  this  type  of  boundary  condition  was  investigated  by  using  various  test  models 
and  a  number  of  useful  conclusions  drawn  through  Monte  Carlo  uncertainty  analysis 

111. 2  Computational  Chemistry:  Inter-atomic  Potential  and  Validation  for  Oxygen-Platinum  Systems 

since  the  most  well-studied  surface  catalysis  problem  in  science  is  the  oxygen-platinum  system, 
the  Molecular  Dynamics  potential  (ReaxFF)  and  simulation  techniques  were  validated 
extensively  with  existing  experimental  and  literature  data 

specifically,  new  simulation  methods  were  developed  using  ReaxFF  to  predict  surface  coverage 
of  oxygen  on  platinum,  as  well  as  model  collisions  of  gas-phase  oxygen  with  platinum  (sticking 
and  scattering  angles),  and  validated  with  high  quality  experimental  data 

111. 3  Quasi-Continuum  Method:  Heat  Transfer  from  Surfaces  into  the  Bulk  Material 

here  the  QC  method  is  extended  to  include  heat  transfer 
unsteady  heat  pulses  are  examined  through  atomistic  systems 

111.4  Computational  Chemistry:  A  New  Finite-Rate  Catalytic  Model  for  Oxygen-Silica 

here  a  finite-rate  catalytic  model  (FRC)  is  developed  from  ReaxFF  MD  simulations  for  oxygen- 
silica  systems 

In  contrast  to  any  prior  work,  all  of  the  following  physical  considerations  are  accounted  for: 

amorphous  bulk  silica  and  resulting  amorphous  surface  structure 
surface  defects  inherent  in  the  amorphous  silica  structure 

defect  structures  and  surface  coverage  resulting  from  exposure  to  high  temperature  oxygen 

catalytic  reaction  mechanisms  specific  to  these  defect  chemistries 

determination  of  activation  energies  and  steric  factors  for  each  catalytic  reaction 

accounting  for  the  effects  due  to  off-normal  collisions  on  reaction  rates 

formulation  of  a  complete  finite-rate  model  ensuring  detailed  balance  for  backwards  rates 

validation  of  the  finite-rate  model  with  existing  experimental  data 

The  result  is  a  new  finite  rate  model  for  air-silica  gas-surface  interactions  that  can  be  used 
directly  with  the  new  CFD  FRC  boundary  condition  developed  in  III.l  above. 
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III.  1  Implementation  of  a  Finite-Rate-Catalytic  Boundary  Condition  in  the  US3D  CFD  Code 

A  finite-rate  catalytic  (FRC)  wall  boundary  condition  has  been  implemented  and  validated  in  the 
US3D  CFD  code  developed  and  maintained  at  the  University  of  Minnesota  [25].  The  user  is  able  to 
specify  a  general  set  of  gas-surface  chemical  reactions  (such  as  those  listed  in  Table  1),  as  well  as  the 
parameters  required  for  the  rates  of  each  reaction  (such  as  those  in  Table  2).  The  implementation  chosen 
for  this  boundary  condition  is  straightforward  but  effective.  The  chemical  rate  equations  are  solved  in  a 
separate  subroutine  and  the  result  is  used  to  simply  re-set  the  species  mass  fractions  in  all  ‘ghost-cells’  in 
a  standard  3D  unstructured  CFD  solver  at  each  time  step.  Thus,  no  modifications  to  the  core  CFD 
solver  itself  are  required,  and  we  believe  this  technique  could  be  adopted  easily  by  other  research 
groups.  We  then  investigated  the  general  behavior  of  this  type  of  boundary  condition  by  using  a  test 
parameterization  for  air-silica  interactions  from  Marschall  et  al.  [5],  resulting  in  a  number  of  important 
conclusions. 


Figure  5  -  Comparison  of  US3D  and  DPLR  results  for  surface  heat  flux  (a)  and  for  surface  coverage  (b). 

Example  results  from  US3D  +  FRC  boundary  condition  are  shown  in  Fig.  5  along  with  results  from 
NASA’s  DPLR  code.  Specifically,  heat  flux  results  for  6  km/s  air- flow  over  a  2m  diameter  cylinder 
(T=200K,  p=0.001kg/m3,  Twan=2250K),  using  the  oxygen-silica  model  from  Table  1  [5]  are  shown  in  Fig. 
5a.  The  new  FRC  wall  boundary  condition  predicts  a  heat  flux  between  the  limiting  assumptions  of  super- 
catalytic  and  non-catalytic  as  expected.  The  corresponding  surface  coverages  (O(s),  N(s),  and  empty  sites 
[s])  around  the  cylinder  are  shown  in  Fig.  5b.  More  importantly,  the  US3D  results  have  been  compared 
with  an  independent  finite-rate  model  implemented  in  the  NASA  DPLR  code.  Both  simulations 
employ  the  same  gas-phase  and  gas-surface  chemical  reactions  and  rates,  and  despite  substantially 
different  numerical  implementations,  the  agreement  is  excellent.  These  results  were  presented  in  a 
series  of  papers  at  the  AIAA  Thermophysics  Conference  in  June  2011  [5,14,26]. 

A  number  of  new  contributions  and  conclusions  resulted  from  this  research,  including: 

1)  Zero-dimensional  gas-surface  simulations  were  used  to  demonstrate  that  the  backwards  rates  of 
surface  recombination  reactions  cannot  be  specified  arbitrarily.  Rather,  if  the  forward  surface 
recombination  rate  is  specified,  the  backwards  rate  is  uniquely  determined  by  the  adsorption/desorption 
equilibrium  constant  as  well  as  the  gas-phase  dissociation  equilibrium  constant.  If  rates  are  not  prescribed 
in  this  manner,  the  system  is  shown  to  drift  from  its  equilibrium  state,  a  situation  that,  physically,  should 
not  be  caused  by  a  catalyst. 

2)  It  was  found  that  under  certain  conditions,  despite  a  constant  surface  temperature,  the  recombination 
efficiency  for  oxygen  atoms  (y0)  increased  by  a  factor  of  4  between  the  stagnation  and  90  degree  points 
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on  the  cylinder.  This  result  can  be  seen  in  Fig.  6a  and  demonstrates  how  recombination  coefficients  may 
not  only  be  dependent  on  surface  temperature,  but  also  on  surface  coverage  which  is  a  non-liner 
result  of  local  gas-phase  concentrations  and  numerous  interacting  gas-surface  mechanisms.  This  is  in 
contrast  to  existing  models  (such  as  shown  previously  in  Fig.  1)  that  model  y  as  a  function  of  temperature 
only. 


(a) 

Figure  6  -  Variation  of  y  due  to  surface  coverage  on  a  constant  temperature  cylinder  (a).  Dominant  reaction 
mechanisms  contributing  to  stagnation  point  heat  flux  determined  from  Monte  Carlo  uncertainty  analysis  (b). 


3)  A  finite-rate  catalytic  model  introduces  a  large  number  of  parameters  that  require  specification,  many 
of  which  are  currently  have  large  uncertainty.  Thus  a  legitimate  concern  involves  uncertainty 
quantification  for  such  a  finite  rate  model  as  a  whole.  To  address  this,  Monte  Carlo  uncertainty 
analysis  was  performed  (involving  >1000  individual  CFD  solutions),  where  each  rate  was  drawn  from  a 
log-normal  distribution  with  roughly  an  order-of-magnitude  variation.  Correlation  of  stagnation  point  heat 
flux  to  each  gas-surface  reaction  rate  was  performed.  Uncertainty  in  the  stagnation  point  chemical  heat 
flux  (<yc)  was  found  to  be  directly  linked  to  uncertainty  in  the  dominant  reaction  rate  at  the  surface 
temperature  considered.  The  progression  of  dominant  reactions  as  wall  temperature  is  increased  is  shown 
in  Fig.  6b,  where  the  reaction  numbers  match  those  in  Table  1.  It  is  our  conclusion  that  careful  analysis 
of  a  specific  gas-surface  interaction  will  often  enable  model  reduction  to  a  few  dominant  reactions 
where  uncertainty  can  be  well  quantified. 

4)  Finally,  at  a  given  surface  temperature,  our  Monte  Carlo  analysis  showed  that  an  increase  in  the 
chemical  heat  flux  (due  to  a  variation  in  the  reaction  rates)  leads  to  a  fraction  of  that  increase  in  the  total 
heat  flux.  This  implies  that  although  an  increase  in  surface  reactivity  increases  the  chemical  heat 
flux,  it  alters  the  boundary  layer  in  a  manner  that  decreases  the  conductive  heat  flux. 
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III. 2  Computational  Chemistry:  Inter-atomic  Potential  and  Validation  for  Oxygen-Platinum  Systems 

The  ReaxFF  Interatomic  Potential  Model  for  Molecular  Dynamics  (MD)  Simulations 

For  our  computational  chemistry  simulations  we  use  the  ReaxFF  potential.  This  potential  is  reactive 
in  nature  and  provides  an  accurate  description  of  dissociation  and  reaction  curves  [27,28].  The  system 
energy  in  ReaxFF  is  divided  into  partial  energy  contributions  similar  to  those  of  empirical  nonreactive 
force  fields  (bond  angle,  torsion,  van  der  Waals,  and  Coulomb)  with  the  distinction  that,  in  ReaxFF,  all 
these  contributions  are  bond  order  dependant  and  the  nonbinded  interactions  (van  der  Waals  and 
Coulomb)  are  included  for  all  atom  pairs  to  ensure  continuity  in  the  energy  description  during  bond 
dissociation: 
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The  bond  order  dependence  allows  a  smooth  transition  from  nonbonded  to  single,  double,  and  triple- 
bonded  systems  by  ensuring  that  the  corresponding  energetic  contributions  disappear  upon  bond 
dissociation.  In  general  ReaxFF  is  trained  to  an  extensive  database  of  DFT  energies,  however,  it  enables 
the  simulation  of  large  atomic  systems  (>1  million  atoms  if  necessary).  Since  the  interatomic  potential 
is  the  fundamental  model  used  in  our  Molecular  Dynamics  (MD)  simulations,  the  ReaxFF 
formulation  and  specific  parameterizations  for  materials  of  interest  must  be  validated  with  theory 
and  experiment  to  the  greatest  extent  possible.  Due  to  the  automotive  industry,  a  large  amount  of 
experimental  data  and  prior  research  focuses  on  oxygen-platinum  catalytic  systems.  Since  platinum 
is  also  used  in  thermocouples  used  to  measure  heat  flux  in  high  enthalpy  facilities,  we  first  validated  our 
numerical  methods  for  oxygen-platinum  systems. 


Molecular  Beam  Experiments  for  Oxygen  Impacting  Pt(l  11) 


The  molecular  dynamics  technique  with  the  ab  initio 
based  classical  reactive  force  field  ReaxFF  was  used 
to  study  the  adsorption  dynamics  of  02  on  Pt(l  1 1)  for 
both  normal  and  oblique  impacts.  Overall,  good 
quantitative  agreement  with  the  experimental  data  was 
found  at  low  incident  energies.  Specifically,  ours  were 
the  first  numerical  simulations  to  reproduce  the 
characteristic  minimum  of  the  trapping  probability  at 
kinetic  incident  energies  around  0.1  eV.  This  feature  is 
well  documented  experimentally  and  our  modeling 
revealed  this  minimum  is  due  to  the  presence  of  a 
physisorption  well  near  the  surface,  combined  with  the 
progressive  suppression  of  a  steering  mechanism 
when  increasing  the  translational  kinetic  energy  or  the 
molecule’s  rotational  energy  because  of  steric 
hindrance.  In  the  energy  range  between  0.1  and  0.4 
eV,  the  sticking  probability  increases,  similar  to 
molecular  beam  sticking  data.  For  very  energetic 
impacts  above  0.4  eV,  ReaxFF  predicted  sticking 
probabilities  lower  than  experimental  sticking  data  by  almost  a  factor  of  3  due  to  an  overall  less  attractive 
ReaxFF  potential  energy  surface  compared  to  experiments  and  density  functional  theory. 


Et  [eV] 

Figure  7  -  Sticking  probability  for  02  collisions 
with  Pt(  111)  in  function  of  impact  energy. 
ReaxFF  result  is  in  black. 


For  oblique  impacts,  the  trapping  probability  was  reduced  by  the  nonzero  parallel  momentum  because  of 
the  PES  corrugation  and  does  not  scale  with  the  total  incident  kinetic  energy.  Furthermore,  our 
simulations  predicted  quasispecular  and  slightly  supraspecular  distributions  of  angles  of  reflection,  in 
accordance  with  molecular  beam  experiments.  Increasing  the  beam  energy  to  1.7  eV  caused  the  angular 
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distributions  to  broaden  and  to  exhibit  a  tail  toward  the  surface  normal  because  molecules  have  enough 
momentum  to  get  very  near  the  surface  and  thus  probe  more  corrugated  repulsive  regions  of  the  PES. 


(a)  (b) 

Figure  8  -  Scattering  angles  for  02  collisions  with  Pt( 111)  at  incident  angles  of  40°  (red)  and  60°  (green).  ReaxFF 
predictions  compared  to  experimental  measurements  for  two  incident  energies:  0.46  eV  (a)  and  1.73  eV  (b). 

Dissociative  Adsorption  of  Oxygen  on  Pt(lll)  and  Surface  Coverage 

In  order  to  investigate  gas-phase  collisions  with  a  material  surface,  one  must  characterize  the  state  of 
the  surface  (the  surface  coverage)  that  exists  under  the  desired  experimental  or  flight  conditions. 
Since  it  is  this  surface  that  gas-phase  molecules  actually  interact  with  under  conditions  of  interest.  We 
developed  a  new  technique  that  combines  the  ReaxFF  potential  with  a  statistical  Monte  Carlo  technique 
to  predict  the  oxygen  coverage  on  a  platinum  surface  in  function  of  gas  temperature  and  pressure  and 
compared  the  predictions  with  extensive  experimental  results. 

The  most  suitable  ensemble  to  study  adsorption  problems  is  the  so-called  grand  canonical  ensemble 
(usually  denoted  with  p,  V,  T),  where  the  temperature  (T),  volume  (V),  and  chemical  potential  (p)  are 
fixed  [29,30].  The  gas  is  assumed  to  be  in  thermal  equilibrium  with  the  surface,  i.e.,  T  =  Ts.  At 
equilibrium,  the  chemical  potential  of  a  gas  X  adsorbed  on  a  surface  M  and  the  chemical  potential  of  X  in 
the  reservoir  are  equal,  specifically: 

Pm-x  =  Fx 

The  Grand  Canonical  Monte  Carlo  (GCMC)  method  is  therefore  an  equilibrium  technique  to  sample 
observables  in  the  p,V,T  ensemble,  for  example  the  average  number  of  particles  adsorbed  on  the  surface 
(or  equivalently  the  surface  coverage  0)  at  a  prescribed  chemical  potential  p. 

The  standard  GCMC  algorithm  [29]  proceeds  through  two  steps: 

(1)  Displacement  of  a  particle:  a  random  particle  at  x  is  selected  and  given  a  random  displacement  5x. 
This  move  is  accepted  with  probability: 

ACC(x  x  +  5x)  =  min(l,  exp{-p[U(x  +  5x)  -  U(x)]}) 

(2)  Insertion  or  removal  of  a  particle:  a  new  particle  is  inserted  at  a  random  position  in  the  volume  V  or  a 
random  particle  in  V  is  removed.  The  insertion  is  accepted  with  probability: 

ACC(N  ->  N+l)  =  min(l,  V/(A3(N  +  1»  *  exp{p[p-U(N+l)+U(N)]})  , 

and  the  removal  is  accepted  with  probability 
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ACC(N  -►  N— 1)  =  min(l,  A3N/V  *  exp{-p[ji+U(N-l)-U(N)]}) . 

Here,  p  =  l/kBT  (kB  is  Boltzmann’s  constant)  and  A  =  [h2/(27imkBT)]°5  is  the  thermal  de  Broglie 
wavelength  (h  is  Planck’s  constant  and  m  is  the  particle  mass).  The  potential  energy  of  the  system  U  is 
obtained  from  the  ReaxFF  potential  and  calculated  by  the  Lammps  MD  package  [31,32]  (distributed  by 
Sandia  National  Lab),  which  can  be  compiled  as  a  set  of  external  routines  called  by  a  driver  program.  The 
use  of  Lammps  allows  the  simulation  of  fairly  massive  systems,  because  the  calculation  of  the  potential 
energy  is  done  in  parallel. 

For  an  ideal  diatomic  gas  X2,  the  chemical  potential  is  written  as: 

P=P°x2  (T,p°)+kBT  ln(px2/p0) 

where  p°x2(T,p0)  is  a  reference  chemical  potential  that  is  only  a  function  of  T,  and  pX2  is  the  pressure  in 
the  reservoir.  p°x2(T,p°)  is  either  tabulated  or  can  be  calculated  with  the  partition  function  of  the  gas  (a 
known  quantity). 
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Figure  9  -  GCMC  convergence  of  surface  coverage  (left).  GCMC  prediction  of  high  oxygen  coverage  and  displaced 
Pt  atoms  displayed  with  experimental  image  (right). 


Thus,  specifying  the  chemical  potential  is  equivalent  to  specifying  the  T,  p  conditions  of  the  gas  next  to 
the  surface.  We  performed  such  GCMC  simulations  for  the  well-studied  problem  of  oxygen  adsorbing  on 
platinum  (111),  for  which  substantial  theoretical,  numerical,  and  experimental  data  exists.  The  results 
were  very  promising,  some  of  which  are  shown  in  Fig.  9.  Specifically,  starting  from  a  vacant  Pt(lll) 
surface,  using  the  GCMC  technique,  equilibrium  surface  coverages  were  obtained  for  a  wide  variety  of  T, 
p  conditions  (even  for  p  approaching  near  vacuum  conditions).  At  these  low  pressures,  the  GCMC 
method  predicted  the  well-known  0.25  mono-layer  (ML)  coverage  and  also  predicted  the  specific  p2x2 
arrangement  of  adsorbed  oxygen  atoms  in  agreement  with  accepted  experimental  evidence.  The  method 
was  also  able  to  predict  high  surface  coverages  where  Pt  atoms  are  actually  displaced  out  of  the  lattice,  in 
qualitative  agreement  with  experimental  images  (seen  in  Fig.  9). 

These  simulations  are  novel  since  they  rely  solely  on  a  reactive  interatomic  potential  (ReaxFF)  to  predict 
the  chemical  structure  of  the  surface  without  pre-specification  of  expected  surface  structures  common  to 
other  methods.  In  addition,  the  technique  is  highly  consistent  in  that  the  same  interatomic  potential  can  be 
used  for  dynamic  MD  simulations.  Therefore  enabling  the  prediction  of  the  in-situ  surface  structure  at 
low  pressures  and  high  temperatures,  followed  by  the  simulation  of  gas  molecules  interacting  with 
this  surface. 
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III. 3  Quasi-Continuum  Method:  Heat  Transfer  from  Surfaces  into  the  Bulk  Material 

In  our  Molecular  Dynamics  simulations  of  gas-surface  interactions  we  thermostat  our  atomic  surface 
system  to  a  certain  temperature  T,  and  assume  that  this  is  representative  of  an  experiment  where  the 
surface  temperature  is  reported  as  the  same  T.  However,  the  question  arises,  is  the  temperature  of  the  first 
few  atomic  layers  (important  for  catalytic  reactions)  actually  the  same  as  the  temperature  of  the  bulk 
material?  More  specifically,  how  is  heat  (perhaps  deposited  due  to  a  catalytic  reaction)  transferred  from 
these  atomic  surface  layers  into  the  bulk?  Secondly,  if  we  can  study  heat  transfer  at  the  atomic  level  as 
phonons  transporting  energy  within  the  solid  material,  then  it  may  be  possible  to  determine  correlations 
between  the  phonon  properties  of  a  material  and  the  vibrational  state  of  molecules  leaving  the  surface. 
This  would  enable  a  much 
more  physical  understanding  of 
the  energy  accommodation 
coefficient  ((3)  that  is  typically 
assumed  to  be  one  for 
hypersonic  flows.  These 
questions  motivated 

fundamental  research  into  heat 
transfer  processes  spanning 
atomistic  to  continuum  scales. 

The  Quasi  Continuum  (QC) 
method  of  co-PI  Tadmor  was 
chosen  as  the  methodology  for 
this  research.  The  QC  method 
rigorously  couples  atomistic 
simulation  to  Finite  Element 
continuum  simulation  for 
multi-scale  material  problems. 

The  QC  method  has  been  applied  to  a  range  of  problems  in  atomic-scale  mechanics,  such  as  fracture  and 
plasticity  at  the  nanoscale.  An  example  of  a  QC  simulation  of  crack  propagation  is  shown  above  (image 
right).  Frame  (a)  shows  the  full  QC  model  corresponding  to  a  domain  of  0.5mmx0.5mm,  most  of  which 
is  treated  by  the  continuum  limit  of  the  method.  Frame  (b)  shows  a  close-up  of  the  tip  of  the  propagating 
crack,  which  is  described  with  full  atomistic  detail.  This  example  shows  the  capability  of  QC  to  correctly 
capture  the  far-field  continuum  behavior  while  still  being  able  to  describe  the  atomistic  processes  of  bond 
breaking  and  dislocation  nucleation  occurring  at  the  crack  tip. 

The  QC  formulation  prior  to  this  grant  was  a  zero-temperature  minimization  scheme.  Formulation  of  a 
finite-temperature,  “hot”-QC  method  has  proven  to  be  a  very  challenging  proposition.  The  key  difficulty 
is  how  to  account  correctly  for  the  entropy  of  the  atoms  that  have  been  integrated  out  of  the  continuum 
region.  Clearly  a  coarse-grained  node  vibrating  in  the  continuum  region  is  not  the  same  thing  as  an  atom 
vibrating  in  the  atomistic  region.  Consequently,  the  naive  methodology  of  simply  “doing  dynamic”  by 
solving  Newton’s  equation  for  the  full  system  (adopted  in  some  multiscale  methods)  cannot  be  correct. 
Our  objective  was  to  use  hot-QC  to  explore  the  effect  of  heat  transfer  between  the  atomistic  region  at  the 
surface  of  the  solid  and  the  bulk  continuum.  This  means  that  hot-QC  must  be  extended  from  a  constant 
temperature  equilibrium  method,  to  one  in  which  a  non-uniform  temperature  field  can  exist. 

In  this  sense  understanding  heat  diffusion  in  the  material  becomes  important  and  the  question  whether 
non-stationary  heat  conduction  is  adequately  modeled  by  Fourier’s  model  or  whether  non-Fourier 
models,  such  as  the  Cattaneo-Vernotte  (CV)  or  Jeffreys  type  (dual-phase-lag)  models,  are  necessary  has 
been  a  matter  of  discussion.  Fourier’s  equation  only  predicts  diffusion  phenomenon  whereas  non-Fourier 


QC  simulation  of  crack  propagation  in  nickel,  (a)  The  full  model.  The  mesh 
has  adapted  to  allow  the  crack  to  grow,  (b)  Detail  of  the  crack  tip.  The  colors 
correspond  to  contours  of  shear  strain  left  in  the  wake  of  dislocations  that  have 
been  emitted  from  the  crack  tip  during  its  motion. 
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models  predict  a  coupled  diffusion-wave  response. 

In  this  work,  using  classical  Non-Equilibrium 
Molecular  Dynamics  (NEMD)  simulations,  we 
investigated  the  process  of  heat  conduction  in  a 
three-dimensional  atomic  beam  of  finite  length.  2] 
Specifically,  we  studied  solid  argon  modeled  using  a  § 
Lennard-Jones  (6-12)  potential  and  silicon  modeled  | 
using  a  Stillinger- Weber  and  Tersoff  potential,  g 
Langevin  and  Nose  Hoover  thermostats  were  used  to 
maintain  a  temperature  gradient  so  that  a  heat  flux 
vector  would  develop  within  the  system.  Results  of 
the  simulations  in  the  form  of  spatio-temporal 
temperature  and  heat  flux  vector  profiles  were 
obtained  and  a  curve  fitting  procedure  based  on  the 
Iteratively  Re-weighted  Least  Squares  (IRLS) 
regression  method  was  used  to  estimate  the 
parameters  appearing  in  the  dual-phase-lag  model. 
The  figure  (image  right)  presents  the  comparison  of 
the  curves  obtained  by  the  continuum  dual-phase-lag 
model  and  unsteady  NEMD  simulations.  It  was 
found  that  the  two  relaxation  times  appearing  in  the 
dual-phase-lag  model  are  almost  equal  and, 
therefore,  Fourier  heat  conduction  is  adequate  for 
modeling  unsteady  heat  conduction  for  insulators  - 
at  least  for  the  cases  studied  here.  Despite  this  result, 
the  use  of  the  dual-phase-lag  model  has  the 
advantage  of  providing  information  on  relaxation 
times  which  cannot  be  obtained  from  the  Fourier 
model  alone.  In  addition,  the  above  procedure  also 
yields  the  thermal  conductivity  of  the  material. 

This  result  was  compared  with  and  found  to  be  in 
agreement  with  two  other  procedures:  Green-Kubo 
and  direct  NEMD.  Of  these  three  approaches,  the 
curve  fitting  procedure  described  above  has  the 
advantage  of  taking  the  least  amount  of  simulation 
time  to  calculate  thermal  conductivity  and 
relaxation  times.  This  is  because  the  procedure  is 
based  on  unsteady  heat  conduction  and  only  a  few 
spatio-temporal  profiles  are  necessary;  in  contrast 
for  the  other  approaches  the  system  needs  to  be 
brought  to  steady  state  which  can  take  a  very  long 
time. 
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Temperature  profiles  obtained  using  NEMD 
simulations  and  a  dual-phase-lag  continuum  model 
for  a  solid  argon  beam.  The  beam  is  initially 
equilibrated  at  10  K  and  the  temperatures  of  the 
ends  are  then  raised  and  held  at  15  K.  The  red 
curves  are  the  averaged  atomistic  results  for  6 
different  initial  conditions  and  the  black  curves  are 
those  obtained  from  dual-phase-lag  model  using 
parameters  obtained  from  the  IRLS  fitting 
procedure.  The  beam  was  divided  into  1 1  bins  and 
Xcm  is  the  position  along  the  beam  of  the  center  of 
mass  of  each  bin  in  Angstroms. 


Temperature  across  length 


The  thermal  conduction  in  amorphous  and  crystalline  silica  material  was  also  studied.  Using  the  NEMD 
direct  method  and  the  curve-fitting  procedure  described  above,  the  thermal  conductivities  of  crystaline 
and  amorphous  silica  at  different  temperatures  have  been  calculated  and  the  results  are  in  agreement 
within  the  range  of  literature  and  experimental  values.  The  figure  (image  above)  shows  the  steady  state 
temperture  profile  for  the  case  where  crystalline  silica  (quartz),  using  the  BKS  potential,  is  maintained  at 
280  K  at  the  ends  of  the  beam  and  320  K  tempertaure  in  the  middle  of  the  beam.  Through  the  curve  fitting 
procedures,  here  too,  it  was  found  that  Fourier  equation  is  sufficient  to  model  heat  transfer. 
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Finally,  studies  have  been  carried  out  to 
understand  how  heat  waves  travel  in  the  atomic 
systems.  The  figure  (image  right)  presents  such  a 
case.  Here,  an  argon  beam  of  size  303  x  4  x  4  is 
divided  into  101  bins,  where  each  bin  consists  of 
243  atoms.  Periodic  boundary  conditions  are 
employed.  A  LJ  (6-12)  pair  potential  is  used  to 
model  the  atomic  interactions.  The  middle  bin  of 
the  beam  is  subjected  to  a  heat  pulse  of  1400 
femtosecond  (fs)  time-duration  with  rise  and  fall 
times  of  200  fs.  Initially,  the  temperature  of  the 
overall  beam  is  OK.  Then,  the  temperature  of  the 
middle  bin  is  increased  to  20K  linearly  during  the 
rise  time  of  200  fs.  After  this,  its  temperature  is 
maintained  at  20K  for  a  duration  of  1000  fs, 
which  is  followed  by  a  linear  downward  decrease 
in  its  temperature  to  OK.  The  rest  of  the  beam  is 

made  to  evolve  under  Newtonian  dynamics  throughout  the  heat  pulse  duration.  The  introduction  of  a  heat 
pulse  generates  several  waves  which  propagate  at  different  speeds  corresponding  to  different  phonon 
modes,  second  sound  waves,  and  diffusive  components.  The  disturbances  are  propagating  in  straight  lines 
and  they  are  labeled  as  LI,  L2,  SI,  HI,  HI,  Hav,  H3,  corresponding  to  different  elastic  and  thermal  waves. 
Ongoing  work  on  the  hot-QC  method  is  repeating  these  caculations  for  crystalline  silica  (quartz)  and 
amorphous  Si02. 
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Heat  wave  disturbances  propagating  through  an 
atomistic  system. 


With  continued  development  of  the  hot-QC  method,  it  should  be  possible  to  investigate  the  heat  pulse 
generated  by  a  surface  reaction,  in  3D  through  an  amorphous  Si02  surface.  Specifically  the  transfer  of 
energy  into  the  bulk  could  be  investigated,  but  also  the  influence  of  such  heat  pulses  on  neighboring 
adsorbed  molecules  would  be  of  interest.  Also,  of  interest  would  be  correlations  between  the  phonon 
characteristics  of  a  given  material  and  the  vibrational  state  of  desorbed  (recombined)  molecules.  If  a 
correlation  exists,  the  implication  would  be  that  the  accommodation  coefficient  ((3)  may  depend  on  the 
material  structure  itself. 
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III.4  Development  of  a  Finite-Rate  Catalytic  Model  for  Oxygen-Silica  using  Computational  Chemistry 

In  this  work  we  now  focus  on  surface  catalysis  of  atomic  oxygen  on  silica  (Si02).  Silica  is  chosen 
because  it  is  a  significant  component  in  both  reusable  (LI900,  LI2200,  FRSI)  and  ablative  (SIRCA) 
thermal  protection  systems  [24].  In  addition,  studies  have  found  that  several  non-Si02  based  thermal 
protection  systems,  such  as  SiC  (at  T  <  1800  K)  and  Ultra  High  Temperature  Ceramics  (ZrB2-SiC  and 
ZrB2-SiC-HfB2  at  T  <  1300  K)  form  thin  Si02  layers  when  exposed  to  air  plasma,  and  act  similarly  to 
pure  silica  from  a  catalytic  perspective  [22,33].  Despite  a  large  body  of  experimental  work  there  is 
uncertainty  in  the  precise  temperature,  pressure,  and  gas  composition  dependence  of  heating  on  silica 
surfaces  due  to  heterogeneous  catalysis  [23].  Reported  values  for  oxygen  recombination  coefficients 
on  silica  in  the  literature  span  several  orders  of  magnitude  (2x1  O'4  <  y  <  4x1  O'1)  at  1000  K  [34-35]. 

The  goal  of  this  work  is  to  develop  a  methodology  for  creating  a  finite  rate  catalytic  model  of  oxygen 
recombination  on  silica  using  molecular  dynamics  simulations  equipped  with  the  ReaxFF  potential.  The 
ReaxFF  potential  is  a  classical  potential  parameterized  from  quantum  chemical  calculations,  which  allows 
for  chemical  reactions  to  occur  over  the  course  of  a  molecular  dynamics  (MD)  simulation.  Using  MD 
simulations,  we  will  identify  the  chemical  structures  on  silica  surfaces  where  recombination  can  occur. 
The  elementary  reactions  in  our  finite  rate  catalytic  model  are  based  on  the  possible  outcomes  of  oxygen 
interaction  with  these  structures.  Molecular  dynamics  simulations  of  individual  events  are  used  to  find  the 
parameters  of  the  functional  forms  of  the  reactions,  including  activation  energies  and  pre-exponential 
factors.  This  methodology  is  used  to  create  a  preliminary  finite  rate  catalytic  model,  which  can  be  used  to 
find  recombination  coefficients  and  surface  coverages  at  a  given  temperature  and  pressure,  and  can  be 
directly  incorporated  into  Computational  Fluid  Dynamics  (CFD)  simulations  using  the  FRC  boundary 
condition  described  in  section  III.  1 . 

Realistic  Amorphous  Si02  Bulk  Structures 

Fundamental  to  any  description  of  recombination  on 
surface  are  the  in-situ  chemical  structures  on  the 
surface  with  which  gas-phase  molecules  interact.  Our 

research  focuses  on  amorphous  silica  (a-Si02)  surfaces, 
has  been  shown  that  Molecular  Dynamics  (MD) 
simulations  with  empirical  potentials  are  able  to  describe 
the  structural  features  on  silica  surfaces  under  standard 
atmospheric  conditions  in  cases  where  comparison  to 
experimental  data  can  be  made  [36,37].  However,  in 
many  experiments  measuring  surface  catalycity,  surfaces 
are  exposed  to  an  air-plasma,  and  there  is  uncertainty  if 
the  same  structures  are  present  under  these  conditions. 

Therefore,  we  will  use  MD  simulations  with  the  ReaxFF 
potential  to  investigate  the  structures  that  occur  on  Si02 
surfaces  before  and  after  they  are  exposed  to  atomic 
oxygen.  All  MD  simulations  are  run  with  the  publicly 
available  molecular  dynamics  program,  LAMMPS, 
which  is  distributed  by  Sandia  National  Labs  [32]. 

Before  modeling  amorphous  silica  surfaces ,  we  had  to  first  validate  that  the  ReaxFF  potential  can 
accurately  reproduce  bulk  amorphous  silica.  To  generate  bulk  a-Si02,  we  performed  molecular 
dynamics  simulations  following  the  annealing  procedure  described  by  Hu  et  al.22  Bulk  (3-cristobalite  was 
given  an  initial  temperature  of  8000  K,  and  propagated  for  20  ps  under  NVT  dynamics  to  randomize  its 
initial  structure.  The  bulk  was  then  cooled  at  50  K/ps  under  NVT  dynamics  until  it  reached  300  K.  The 
system  was  propagated  for  a  further  40  ps  under  NPT  dynamics  at  300  K,  1  atm,  with  the  last  20  ps  used 
to  collect  statistics  about  the  structure.  We  used  11,616  atoms  in  our  simulations,  enough  to  ensure 


Figure  10  -  Amorphous  Si02  bulk  structure 
created  by  annealing  (3-cristobalite  using 
ReaxFF. 
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sufficient  statistics  about  the  structure  were  collected  [38].  Simulations  were  performed  with  the  ReaxFF 
potential  and  the  BKS  potential  using  the  modifications  outlined  in  the  work  by  Jee  et  al.  [39].  An  image 
of  the  resulting  structure  is  seen  in  Fig.  10. 

The  total  correlation  function  T(r)  can  be  used  to 
directly  compare  the  computationally  generated 
annealed  bulk  structure  to  neutron  scattering 
experiments.  This  function  is  generated  from  the 
partial  radial  distribution  functions  as  described 
by  Nakano  et  al.[40]  As  shown  in  Fig.  11,  the 
total  correlation  function  of  annealed  a-Si02 
created  with  both  potentials  is  in  relatively  good 
agreement  with  experimental  results  [41]. 

However,  the  ReaxFF  potential  under-predicts  the 
location  and  magnitude  of  the  secondary  peak, 
which  represents  the  average  distance  between  O- 
O  nearest  neighbors,  as  well  as  several  tertiary 
peaks.  The  BKS  potential  is  in  better  agreement 
with  experimental  measurements  on  bulk  a-Si02, 
however,  this  potential  lacks  the  ability  to 
describe  chemical  reactions  such  as  gas  phase 
oxygen-oxygen  bonding  because  of  its  two-body 
nature.  Ref.  [42]  contains  detailed  comparisons  of 
ReaxFF  predictions  to  an  entire  set  of 
experimental  measurements  including 
coordination  numbers  and  angle  distributions.  The  ReaxFF  potential  was  found  to  be  in  quite  good 
agreement  with  experimental  measurements,  which  tended  confidence  that  the  ReaxFF  potential 
could  adequately  model  bulk  a-Si02;  a  pre-requisite  for  modeling  a-Si02  surfaces. 

Realistic  Amorphous  Si02  Surface  Structures 

The  simplest  way  to  generate  an  a-Si02  surface  is  to  cleave  the  previously  generated  bulk  along  a 
convenient  plane.  However  this  method  creates  an  idealized  surface  that  is  terminated  by  a  number 
of  highly  reactive  broken  bonds.  To  generate  more  realistic  surfaces,  we  performed  MD  annealing 
simulations  on  an  a-Si02  surface  cleaved  from  the  previously  generated  bulk.  To  create  surfaces,  we 
followed  the  procedure  detailed  in  Fogarty  et  al.  [43]  A  slab  of  cleaved  a-Si02  was  heated  at  25K/ps  from 
300  K  to  4000  K  under  NVT  dynamics  and  then  cooled  back  to  300  K  at  same  rate.  The  system  was 
propagated  for  an  additional  40  ps  under  NPT  dynamics  at  300  K,  1  atm  with  the  final  20ps  used  to 
collect  statistics.  All  slabs  were  20  Angstroms  thick  with  areas  chosen  so  that  surfaces  had  a  statistically 
significant  number  of  defects,  as  described  in  the  upcoming  section.  In  general,  we  found  that  longer 
annealing  simulations,  or  additional  annealing  cycles,  produced  surfaces  with  slightly  lower  surface 
energies  and  slightly  fewer  defects.  However,  we  found  that  the  same  types  of  defects  occur  on  surfaces 
annealed  for  longer  times.  Therefore  we  are  confident  that  although  surfaces  may  contain  more  defects 
that  they  would  if  annealed  at  slower  and  more  experimentally  realistic  rates,  they  still  contain  the  salient 
structural  features  of  an  a-Si02  surface.  We  then  used  the  annealed  a-Si02  surfaces  as  a  basis  for  realistic 
amorphous  silica  surfaces  that  we  will  subsequently  expose  to  atomic  oxygen. 

Catalytic  Defects  on  Amorphous  Si02  Surfaces  After  Exposure  to  Dissociated  Oxygen 

Unlike  surfaces  created  from  crystalline  bulk  materials,  annealed  amorphous  Si02  surfaces  contain 
a  finite  number  of  catalytic  defects  (surface  sites  where  catalytic  events  may  be  energetically 
favorable).  However,  these  defects  may  change  in  number  or  in  chemical  structure  when  exposed  to  a 


Figure  1 1  -  Total  correlation  function  T(r)  plotted  for 
both  ReaxFF  and  BKS  potentials  compared  with 
experiment  [41]. 
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high  temperature  dissociated  gas,  as  may  be  the  case  in  air  plasma  experiments  and  hypersonic  flight.  In 
order  to  simulate  exposure  to  atomic  oxygen,  we  used  MD  simulations  with  a  novel  Flux  Boundary 
Condition  (FBC);  a  method  that  exposes  a  surface  to  an  ideal  gas  at  a  given  temperature,  pressure, 
and  composition.  A  more  detailed  description  of  this  method  can  be  found  in  Ref.  [42].  In  Flux 
Boundary  Condition  simulations,  atoms  are  generated  at  random  points  on  a  plane  above  the  surface  with 
a  frequency  corresponding  to  the  flux  of  an  ideal  gas  through  that  plane.  This  plane  is  10  Angstroms 
above  the  surface,  beyond  the  force  cutoff  of  the  potential.  Impinging  atoms  are  given  velocities  sampled 
from  a  Maxwell-Boltzmann  distribution  and  a  random  incident  angle.  Atoms  or  molecules  moving  away 
from  the  surface  at  a  distance  greater  than  10  Angstroms  are  deleted,  allowing  us  to  simulate  only  the  gas- 
surface  interface  region,  and  gas  phase  collisions  in  this  region  are  rare.  A  diagram  of  this  method  is 
shown  in  Fig  12.  Over  the  course  of  a  simulation,  gas  atoms  adsorb  on  the  initially  empty  surface,  which 
eventually  reaches  a  steady  state  population  of  adsorbed  atoms.  At  steady  state  we  can  determine  the 
concentration  and  nature  of  chemical  defects  for  an  amorphous  Si02  surface  exposed  to  a 
dissociated  gas  at  a  given  temperature  and  pressure. 
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Figure  12  -  Schematic  of  the  Flux  Boundary  Condition  simulation  technique  (red=oxygen,  yellow=silicon) 

We  are  primarily  interested  in  structures  on  the  surface  where  recombination  can  occur.  The  most 
prevalent  structures  are  bridging  oxygen  atoms  (Si  -  O  -  Si) .  These  bridging  oxygen  atoms  are  typical  of 
the  stable  reconstructions  shown  much  earlier  in  Fig.  3a,  and  are  strongly  bound  to  the  surface  (-10  eV 
bonding  energy).  The  energy  released  in  gas  phase  oxygen  recombination  is  only  5.16  eV,  so  the  direct 
recombination  of  an  impinging  gas-phase  oxygen  atom  with  a  bridging  oxygen  atom  is  very  unlikely,  and 
thus,  these  stable  reconstructions  are  non-catalytic.  Therefore,  we  expect  that  recombination  will 
occur  at  defects  on  the  surface.  We  find  that  the  most  populous  defects  on  surfaces  exposed  to  atomic 
oxygen  are  the  non-bridging  oxygen  atom  (=Si-0*)  (Fig.  13b)  and  adsorbed  molecular  oxygen 
(=  Si  -  02)  (Fig.  13c).  We  also  find  that  the  annealed  a-Si02  surfaces  have  a  significant  number  of  under¬ 
coordinated  silicon  atoms  (=Si*)  (Fig.  13a)  before  exposure  to  the  dissociated  gas,  and  also  during 
exposure  when  adsorbed  oxygen  atoms  are  removed  due  to  catalytic  processes. 

Visualizations  of  Si02  surfaces  with  these  defects  highlighted  are  shown  in  Fig.  14.  This  is 
representative  of  the  actual  surface  that  gas-phase  atoms  interact  with.  Figure  14  corresponds  to  a  high 
pressure,  however,  we  have  performed  numerous  studies  showing  that  at  lower  pressures,  the 
concentration  of  defects  is  reduced,  but  the  defect  chemistries  themselves  remain  identical.  There  is 
experimental  evidence  for  the  existence  of  the  (=Si-)  and  (=Si-0  )  defects  on  silica  surfaces 
under  irradiation  and  fracture.  The  (=Si*)  defect  has  been  extensively  studied,  and  has  been  identified 
on  irradiated  and  vacuum  fractured  quartz  and  amorphous  silica  using  Electron  Spin  Resonance  (ESR) 
[44,45].  ESR  in  conjunction  with  isotope  effects  have  been  used  to  identify  the  (=Si-0*)  defect  on 
irradiated  amorphous  Si02.  The  (=  Si  * )  and  (=  Si  -  O  • )  structures  have  also  been  observed  in  other  MD 
simulations  of  a-Si02  surfaces  simulated  using  different  interatomic  potentials  [36,46,47].  We  are  not 
aware  of  any  experimental  evidence  of  the  (=Si-02)  structure,  however  recent  DFT  calculations 
performed  by  the  Truhlar  group  at  the  University  of  Minnesota  as  part  of  a  new  AFOSR  MURI  project 
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confirm  the  energetics  and  stability  of  this  defect.  To  study  recombination  on  realistic  silica  surfaces,  our 
finite  rate  catalytic  model  thus  focuses  on  the  (=  Si  * ) ,  (=  Si  -  O  • ) ,  and  (=  Si  -  02)  defects  (Fig.  13). 


(a)  (-  Si-)  (b)  (  =  Si-0-)  (c)  OSi-02) 

Figure  13  -  Defects  observed  on  annealed  amorphous  Si02  surfaces  exposed  to  dissociated  oxygen. 


Figure  14  -  Visualization  of  annealed  amorphous  Si02  after  exposure  to  dissociated  oxygen  at  high 
pressure.  Defects  are  highlighted.  Blue  =  (=  Si  -  O  * ) .  Green  =  (=  Si  -  02) . 


Gas-Phase  Collisions  with  Defect  Surfaces 
Including  Off-Normal  Collisions 

For  our  finite  rate  catalytic  model,  we  considered 
direct  gas  phase  interactions  with  the  (=Si  ), 
(=Si-0  ),  and  (=Si-02)  defects.  The  rates  of 
reactions  at  these  defects  were  found  with  MD 
simulations  of  the  collision  of  single  gas-phase 
atomic  oxygen  atoms  with  the  defects.  Simulations 
were  performed  for  isolated  defects  on  an  otherwise 
reconstructed  surface.  It  can  be  shown  that  on  the 
scale  of  MD  simulations  (Angstroms)  the  defect 
density  is  low  and  gas  phase  collisions  occur  rarely 
and  are  spread  over  large  surface  areas.  Thus  defects 
do  not  interact  with  other  defects  and  sequential  gas- 
phase  collisions  do  not  interact  with  one  another.  The 
relative  energies  of  defect  structures  on  the 
reconstructed  surface  are  shown  in  Fig  15. 


Figure  15  -  Relative  energy  of  defects. 
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A  set  of  reactions  was  chosen  based  on  the  range  of  possibilities  for  gas  phase  oxygen  atoms 
interacting  with  defects  on  the  surface,  as  shown  in  Table  3.  The  first  reaction  is  atomic  oxygen 
adsorption,  which  is  considered  irreversible  based  on  the  high  energy  of  adsorption  (-4.21  eV).  We  have 
found  that  oxygen  atoms  colliding  with  the  (=Si-0*)  site  can  either  form  (=Si-02)  (reaction  2),  or 
recombine  via  Eley-Rideal  (ER)  recombination  (reaction  3).  Additionally,  we  have  observed  that  atomic 
oxygen  atoms  encountering  (=Si-02)  can  either  replace  the  02  molecule  or  combine  with  one  of  the 
oxygen  atoms  in  the  (=  Si  -  02)  structure  via  ER  recombination  (reaction  4).  In  either  case  the  end  result 
is  identical:  an  absorbed  atomic  oxygen  and  a  gas  phase  molecular  oxygen.  We  also  consider  the 
possibility  of  desorption  of  molecular  oxygen  from  the  (=  Si  -  02)  site  (reaction  5).  The  variables  used  in 
the  finite  rate  catalytic  model  to  describe  the  concentrations  of  defects  and  surface  coverages  are  shown  in 
Table  4.  The  total  number  of  active  sites  [S]  is  constant  and  is  equal  to  the  sum  of  the  concentrations  of 
individual  defects. 


Mechanism 

Name 

0  +  =Si-  -►  =Si-0- 

Atomic  Adsorption 

(i) 

0  +  =Si-0-  — >  =Si-02 

02  Formation 

(2) 

0  +  =Si-0-  -►  02  +  =Si- 

ER  Recombination 

(3) 

0  +  =Si-02  — ►  02  +  =Si-0' 

ER  Recombination  11/  02  Replacement 

(4) 

=Si-02  — >  02  +  =Si- 

02  Desorption 

(5) 

Table  3  -  Elementary  gas- surface  interactions  forming  the  finite-rate  catalytic  model 


Surface  Concentrations 

[0]8  Concentration  of  adsorbed  oxygen  (=Si-0  )  on  the  surface  (m~2) 

[02]s  Concentration  of  adsorbed  molecular  oxygen  (=Si-02)  on  the  surface  (m-2) 
[E]s  Concentration  of  vacant  sites  (=Si*)  sites  on  the  surface  (m-2) 

[S]  Total  concentration  of  active  sites  on  the  surface  (m-2) 

(by  definition  [0\s  +  [02}s  +  [E\s  =  [5]) 

Fractional  Coverages 
Go  [O],  /  [5] 

9o,  [02],  /  [5] 

Oe  [E\.  /  [5] 

(by  definition  do  +  Qo$  +  Or  =  1) 

Table  4  -  Surface  coverage  variables  completing  the  finite-rate  catalytic  model 
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Adsorption:  (O  +  =Si*  =>  =  Si  -  O  •) 

The  energy  of  adsorption  of  oxygen  is  -4.21  eV,  as  shown  in  the  potential  energy  surface  (PES)  in  Fig 
16a.  Given  this  deep  potential  well,  we  assume  that  the  sticking  probability  for  O  on  this  site  is  one. 


(a)  PES  for  atomic  oxygen  adsorption  on  =Si- 


' V 


(b)  Diagram  of  atomic  oxygen  adsorption  on  =Si- 


Figure  16  -  Atomic  oxygen  adsorption  on  (=  Si  • ) 


Molecular  Oxygen  Formation  and  ER-Recombination: 

(O  +  =Si-0 •  =>  =  Si  -  02)  and  (O  +  =  Si  -  O •  =>  02  +  =  Si*) 


0.16 


(a)  PES  for  O2  formation  (b)  Diagram  of  O2  formation  simulations 


Figure  17  -  Molecular  oxygen  formation  (O  +  =  Si  -  O  *  =>  =  Si  -  02) 


As  shown  in  the  PES  in  Fig.  17a,  there  is  a  well-defined  transition  state  for  02  formation.  To  find  the  pre¬ 
exponential  factor  and  activation  energy  for  this  reaction,  we  performed  molecular  dynamics  simulations. 
An  oxygen  atom  was  placed  10  Angstroms  above  the  surface  and  given  a  velocity  sampled  from  the 
Maxwell-Boltzmann  distribution.  The  velocity  of  the  atom  was  normal  to  the  surface  (off-normal 
collisions  will  be  addressed  in  an  upcoming  section  of  the  report).  The  atom  was  placed  directly  above  the 
(=Si-0*)  defect  and  then  given  a  random  displacement  parallel  to  the  surface  between  0  and  2.5 
Angstroms  to  account  for  the  radius  of  the  site.  To  keep  the  surface  from  translating  over  the  course  of 
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MD  simulations,  the  bottom  atomic  layer  was  frozen.  The  two  atomic  layers  above  the  frozen  layer  were 
thermalized  with  the  Langevin  thermostat  in  order  to  maintain  the  temperature  of  the  surface.  The  top 
seven  atomic  layers  were  allowed  to  move  freely  so  as  not  to  influence  the  gas-surface  chemistry. 
Simulations  were  propagated  until  the  atom  collided  and  possibly  reacted  with  the  defect  on  surface.  Over 
2000  trajectories  were  carried  out  at  each  temperature,  with  temperatures  ranging  from  250-1750  K  with 
an  interval  of  250  K.  The  results  of  these  simulations  are  shown  in  Fig.  18a,  where  error  bars  are  based  on 
the  Wald  method  with  a  95%  confidence  interval.  At  all  temperatures,  the  most  probable  reaction  is  02 
formation.  A  log  plot  of  the  probability  of  02  formation,  as  seen  in  Fig.  18b,  shows  that  the  probability  of 
this  reaction  follows  an  exponential  trend,  with  activation  energy  Ea  =  0.149  eV,  in  agreement  with  the 
PES  in  Fig.  17a.  From  the  log  plot,  we  also  find  that  the  pre-exponential  factor  for  this  reaction  is  0.85. 


Figure  18  -  Probability  of  02  formation  and  ER-Recombionation  (a).  Activation  energy  for  02  formation  (b). 


ER-Recombination/O 2  Replacement  { O  +  =Si  02=>  02  +  =Si-0*) 


Simulations  for  this  reaction  were  carried  out  in  an  identical  manner  to  those  for  02  formation  as  shown  in 
Fig.  19.  Over  4000  trajectories  are  carried  out  at  each  temperature.  Although  there  are  two  distinct 
processes,  for  the  purposes  of  the  finite  rate  catalytic  model  we  fit  the  sum  of  both  processes  because  they 
produce  the  same  result.  As  shown  in  Fig.  19b  this  results  in  an  activation  energy  of  Ea  =0.457  eV  and  a 
pre-exponential  factor  of  0.717. 


Figure  19  -  02  Replacement  simulation  (a).  Activation  energy  for  02  Replacement  +  ER-Recombination  (b). 
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Desorption  of  02  (=Si  *  02  =>  02  +  =Si*) 

Our  results  indicate  that  there  is  no  significant  saddle  point  in  the  PES  for  02  desorption.  According  to  the 
ReaxFF  potential,  this  PES  does  not  significantly  change  for  rotational  conformations  of  the  02  molecule 
in  plane  with  the  surface.  Direct  MD  simulations  of  the  desorption  of  02  molecules  from  this  site  would 
take  computationally  inaccessible  times  due  to  the  high  energy  of  desorption.  According  to  transition  state 
theory,  the  pre-exponential  factor  for  the  desorption  of  an  immobile  species  through  an  immobile 
transition  state  is  1012  -  4xl013  [48].  Therefore  we  use  a  pre-exponential  factor  of  1013  and  an 
activation  energy  Ea=  2.88  eV  (as  shown  in  Fig.  15)  as  parameters  for  the  finite  rate  catalytic  model. 

Accounting  for  Off-Normal  Collisions 

We  then  analyzed  the  effect  of  off-normal  collision  angles  on  the  specific  reactions.  In  previous 
simulations,  atoms  were  placed  above  the  surface  randomly  within  a  disk  with  a  radius  of  2.5  Angstroms. 
For  off-normal  collisions,  we  applied  a  rotation  matrix  to  the  disk  and  the  initial  velocity  vector  of  the 
impinging  atom.  This  is  depicted  in  Fig.  20  for  the  02  formation  reaction. 

Figure  20  shows  the  trajectories  of  many  single  collision  simulations.  In  Fig.  20a,  at  0=30°,  impinging 
atoms  maintain  their  disk-like  spatial  distribution  until  the  impact  at  the  site,  where  they  either  react  or 
scatter.  In  Fig.  20b,  at  0=75°,  there  is  a  slight  distortion  of  the  disk  as  it  nears  the  surface  because  the 
surface  is  repulsive.  By  increasing  the  angle  beyond  70°,  we  increase  the  number  of  atoms  that  enter  a 
repulsive  area  above  the  surface  before  hitting  the  defect,  decreasing  the  number  of  reacting  atoms  due  to 
the  fact  that  atoms  placed  in  the  repulsive  region  almost  instantly  scatter.  In  reality,  atoms  approaching  at 
such  angles  would  have  encountered  the  repulsive  region  and  scattered  earlier  in  their  trajectories. 


(a)  <£  =  30  (b)  <p-15 

Figure  20  -  Images  from  off-normal  collision  simulations  of  the  02  formation  reaction. 

Figure  21  clearly  demonstrates  that  the  angle  does  affect  the  outcome  of  the  02  formation  and  02 
replacement  reactions.  Both  these  reactions  are  dominant  mechanisms  for  recombination.  Thus  we 
conclude  that  off-normal  collisions  must  be  accounted  for  when  parameterizing  recombination 
reaction  rates  based  on  molecular  dynamics  simulations. 
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Figure  21  -  Reaction  probabilities  in  function  of  impact  angle. 


The  Complete  Finite  Rate  Catalytic  Model  Predicted  by  ReaxFF 

The  complete  model,  which  also  includes  low  temperature  LH  reactions,  is  completely  detailed 
below.  When  used  to  determine  recombination  coefficients  in  function  of  temperature,  the  new  FRC 
model  can  be  compared  with  available  experimental  data  as  seen  below  in  Fig.  22.  Unlike  previous 
researchers  claims,  we  strongly  argue  that  no  model  developed  from  computational  chemistry  is  able  to 
precisely  predict  the  magnitude  of  experimentally  inferred  y  values.  Since  these  values  should  scale 
directly  with  the  number  of  active  catalytic  sites  (the  defect  density)  at  the  macroscale ,  which  may  be 
highly  dependent  on  the  material  quality,  surface  roughness,  and  the  gas-phase  diffusion  model  used  to 
interpret  the  experimental  observations  (varies  for  different  experiments).  None  of  these  real  physical 
considerations  can  be  taken  into  account  in  atomistic  MD  simulations.  However,  we  argue  that  the  trend 
of  y  with  temperature  should  be  predicted  from  computational  chemistry.  This  trend  is  completely 
dependent  on  the  dominant  chemical  reaction  mechanisms,  their  pre-exponential  factors,  and  activation 
energies;  all  of  which  are  predicted  by  our  simulations.  In  Fig.  22,  the  density  of  active  surface  sites  in  our 
FRC  model  was  set  to  match  the  magnitude  of  y  found  by  experiment  [C].  Although  experimental 
magnitudes  do  not  agree  with  each  other,  at  high  temperatures  the  slope  of  y  vs.  T  is  somewhat  consistent 
and  our  FRC  model  lies  within  the  experimental  uncertainty. 


O  +  Ec  v=^  Oc 

Atomic  Oxygen  Adsorption 

(i) 

O  +  Oc  02  +  Ec 

Eley  Rideal  Recombination 

(2) 

O  +  Oc  02c 

Molecular  Oxygen  Formation 

(3) 

O  +  02c  02  -I-  Oc 

Molecular  Oxygen  Replacement 

(4) 

02c  ^  02  +  Ec 

Molecular  Oxygen  Desorption 

(5) 

O  +  Ef  O  f 

Atomic  Oxygen  Physisorption 

(6) 

Of  +  Ec  E f  +  Oc 

Physisorbed  to  Chemisorbed 

(7) 

Of  +  Oc  E f  +  Ec  +  02 

Langmuir  Hinshelwood  Recombination 

(8) 

[O]  Gas  Phase  Atomic  Oxygen  Concentration  m-3 

[02]  Gas  Phase  Molecular  Oxygen  Concentration  m-3 

[Ec]  Empty  Chemisorption  Site  Density  m-2 

[E /]  Empty  Physisorption  Site  Density  m-2 

[Ocj  Chemisorbed  Atomic  Oxygen  Density  m-2 

[02c]  Chemisorbed  Molecular  Oxygen  Density  m-2 

[O /]  Physisorbed  Atomic  Oxygen  Density  m-2 

rs  Site  Radius  (m) 

cx  Average  speed  (m/s) 
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Figure  22  -  Comparison  of  the  recombination  coefficients  predicted  by  the  FRC  model  compared  to  experiment. 
Experimental  data  is  taken  directly  from  the  published  articles  listed  below. 

[A]  Balat-Pichelin,  M.,  Badie,  J.,  Berjoan,  R.,  and  Boubert,  P.,  “Recombination  coefficient  of  atomic 
oxygen  on  ceramic  materials  under  earth  re-entry  conditions  by  optical  emission  spectroscopy,”  Chemical 
Physics,  Vol.  291,  No.  2,  2003,  pp.  181-194. 

[B]  Stewart,  D.  A.,  “Surface  Catalysis  and  Characterization  of  Proposed  Candidate  TPS  for  Access-to- 
Space  Vehicles,”  NASA  Technical  Memorandum  1 12206,  1997. 

[C]  Kim,  Y.  C.  and  Boudart,  M.,  “Recombination  of  O,  N,  and  H  Atoms  on  Silica:Kinetics  and 
Mechanism,”  Langmuir,  Vol.  7,  1991. 

[D]  Marschall,  J.,  “Experimental  Determination  of  Oxygen  and  Nitrogen  Recombination  Coefficients  at 
Elevated  Temperature  Using  Laser-Induced  Fluorescence,”  AIAA,  Baltimore,  MD,  1997. 
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